Problème torsion d’arbres#

Librairies#

# Importation des librairies
import dolfinx  # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type  # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem  # Résolution de systèmes linéaires
from mpi4py import MPI  # Interface pour le calcul parallèle
import ufl  # Unified Form Language pour les formulations variationnelles
import numpy as np  # Calcul numérique efficace
import matplotlib.pyplot as plt
import math  # Module de fonctions mathématiques standard de Python
import pyvista as pv  # Visualisation 3D scientifique
import gmsh  # Générateur de maillage 3D
import meshio  # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile  # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem  # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh  # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc  # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
Hide code cell output

Géométrie#

Définition de la géométrie et du maillage#

import gmsh
import numpy as np
import meshio
from dolfinx.io import XDMFFile
from mpi4py import MPI
import ipywidgets as widgets
from IPython.display import display
import pyvista as pv
from dolfinx import plot

# Initialisation du plotteur PyVista en dehors de la fonction
p = pv.Plotter(notebook=True)
p.set_background("grey")

# Fonction pour créer et visualiser la géométrie
def create_ellipse(rho_x, rho_y, h=1.0, lc=0.02):
    global domain
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add("ellipse")

    # Création de l'ellipse
    ellipse = gmsh.model.occ.addEllipse(0, 0, 0, rho_x, rho_y)
    curve_loop = gmsh.model.occ.addCurveLoop([ellipse])
    surface = gmsh.model.occ.addPlaneSurface([curve_loop])

    # Extrusion pour obtenir un cylindre elliptique
    gmsh.model.occ.extrude([(2, surface)], 0, 0, h)
    gmsh.model.occ.synchronize()

    # Configuration du maillage
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
    gmsh.model.mesh.generate(3)

    # Sauvegarde et conversion du maillage
    gmsh.write("ellipse.msh")
    msh = meshio.read("ellipse.msh")
    meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))
    gmsh.finalize()

    # Lecture du maillage avec FEniCSx
    with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
        domain = xdmf_file.read_mesh(name="Grid")
        domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
    
    # Conversion du maillage pour la visualisation avec PyVista
    u_topology, u_cell_types, u_geometry = plot.vtk_mesh(domain)
    u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)

    # Effacer la scène précédente et ajouter le nouveau maillage
    p.clear()
    # Ajout du maillage à la scène de visualisation
    p.add_mesh(u_grid, 
               show_edges=True,  # Affiche les arêtes du maillage
               scalar_bar_args={
                   "title": "u",  # Titre de la barre de couleur
                   "title_font_size": 24,
                   "label_font_size": 22,
                   "shadow": True,
                   "italic": True,
                   "font_family": "arial",
                   "vertical": False  # Orientation horizontale de la barre de couleur
               })

    # Ajout d'un titre à la visualisation
    p.add_text("Cylindre Mesh", font_size=24, color="black", position="upper_edge")

    # Ajout des limites de la boîte englobante
    p.show_bounds(color="black")

    # Ajout des axes de coordonnées
    p.add_axes(color="black")

    # Définition de la couleur de fond
    p.set_background("grey")
    
    # Mise à jour de la scène
    p.show()
    
    return rho_x, rho_y, h
    
# Widgets interactifs pour a, b, et h
rho_x_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-grand axe')
rho_y_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-petit axe')
h_slider = widgets.FloatSlider(value=1.0, min=0.5, max=3.0, step=0.1, description='Hauteur')

# Interface utilisateur et fonction interactive
ui = widgets.VBox([rho_x_slider, rho_y_slider, h_slider])
out = widgets.interactive_output(create_ellipse, {'rho_x': rho_x_slider, 'rho_y': rho_y_slider, 'h': h_slider})

display(ui, out)

# Fonction pour récupérer les valeurs actuelles des curseurs
def get_slider_values():
    current_rho_x = rho_x_slider.value
    current_rho_y = rho_y_slider.value
    current_h = h_slider.value
    return current_rho_x, current_rho_y, current_h

# Exemple d'utilisation : récupération des valeurs
rho_x, rho_y, h = get_slider_values()
rho = np.array([rho_x, rho_y])  # Demi-axes dans les directions x et y

Configuration du problème#

Définition de l’espace fonctionnel#

# Définition de l'espace fonctionnel pour le problème
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

Définition des frontière du domaine#

Hide code cell content
# Définition des fonctions pour identifier les différentes faces du cylindre
def down_face(x):
    """Identifie la face inférieure du cylindre (z = 0)"""
    return np.isclose(x[2], 0)

def top_face(x):
    """Identifie la face supérieure du cylindre (z = h)"""
    return np.isclose(x[2], h)

def lateral_face(x):
    """
    Identifie la surface latérale du cylindre elliptique
    Utilise l'équation de l'ellipse : x^2/a^2 + y^2/b^2 = 1
    """
    tolerance = 1e-5  # Tolérance pour la comparaison numérique
    return (np.isclose((x[0]**2 / rho[0]**2 + x[1]**2 / rho[1]**2), 1.0, atol=tolerance)) & (0 <= x[2]) & (x[2] <= h)

# Dimension des facettes (2D pour un domaine 3D)
fdim = domain.topology.dim - 1

# Localisation des entités (facettes) correspondant à chaque face
Sigma_l = mesh.locate_entities(domain, fdim, lateral_face)
Sigma_0 = mesh.locate_entities(domain, fdim, down_face)
Sigma_h = mesh.locate_entities(domain, fdim, top_face)

# Préparation des marqueurs de facettes pour les conditions aux limites
marked_facets = np.hstack([Sigma_l, Sigma_0, Sigma_h])
marked_values = np.hstack([
    np.full_like(Sigma_l, 1),  # Marqueur 1 pour la face latérale
    np.full_like(Sigma_0, 2),  # Marqueur 2 pour la face inférieure
    np.full_like(Sigma_h, 3)   # Marqueur 3 pour la face supérieure
])

# Tri et création des tags de maillage
sorted_indices = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, 
                                  marked_facets[sorted_indices], 
                                  marked_values[sorted_indices])

# Définition de la mesure pour les facettes marquées
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)


# Création d'un dictionnaire pour mapper les noms aux valeurs numériques
surface_ids = {'Sl': 1, 'S0': 2, 'Sh': 3}

Visualisation des frontière du domaine#

Hide code cell outputs
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Conditions aux limites#

# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)

# Application des conditions aux limites sur chaque face
bc0_S0 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V)  # Face inférieure
bc0_Sh = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V)  # Face supérieure
bc0_Sl = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V)  # Surface latérale

Propriétés matériau#

import ipywidgets as widgets
from IPython.display import display, clear_output
from ipywidgets import HBox, VBox

# Widgets pour les constantes élastiques en GPa
mu_input = widgets.FloatText(value=80, description="μ = ", step=1)
mu_label = widgets.Label(value="GPa")

lambda_input = widgets.FloatText(value=120, description="λ = ", step=1)
lambda_label = widgets.Label(value="GPa")

# Organisation avec les labels d'unité
mu_box = HBox([mu_input, mu_label])
lambda_box = HBox([lambda_input, lambda_label])

# Fonction pour obtenir les constantes élastiques en Pa
def display_elastic_constants(change):
    mu = mu_input.value * 1e9  # Conversion en Pa
    lambda_ = lambda_input.value * 1e9  # Conversion en Pa
    
    # Afficher les valeurs en Pa
    clear_output(wait=True)  # Clear previous output
    
    return mu, lambda_  # Retourner les valeurs en Pa

# Liaison de la fonction d'affichage aux événements de changement de valeur des widgets
mu_input.observe(display_elastic_constants, names='value')
lambda_input.observe(display_elastic_constants, names='value')

# Affichage des widgets
ui_elastic = VBox([mu_box, lambda_box])
display(ui_elastic)
mu, lambda_ = display_elastic_constants('value')
print(f"Module de cisaillement (μ) en Pa: {mu:.2f}")
print(f"Premier paramètre de Lamé (λ) en Pa: {lambda_:.2f}")
Module de cisaillement (μ) en Pa: 80000000000.00
Premier paramètre de Lamé (λ) en Pa: 120000000000.00

Définition des chargements#

import ipywidgets as widgets
from IPython.display import display
import numpy as np

def create_angle_widget():
    # Création du widget curseur
    angle_slider = widgets.FloatSlider(
        value=20.0,  # Valeur initiale à 20 degrés
        min=0,
        max=360,
        step=1,
        description='Angle (°):',
        disabled=False,
        continuous_update=True,
        orientation='horizontal',
        readout=True,
        readout_format='.0f',
    )

    # Widget de sortie pour afficher les valeurs
    output = widgets.Output()

    # Fonction pour mettre à jour l'affichage
    def update_angle(change):
        ang_degrees = change['new']
        ang_radians = np.radians(ang_degrees)
        
        with output:
            output.clear_output(wait=True)
            print(f"Angle : {ang_degrees:.0f}° ({ang_radians:.2f} radians)")
        
        return ang_radians

    # Attacher la fonction de mise à jour au widget
    angle_slider.observe(update_angle, names='value')

    # Afficher le widget et la sortie
    display(widgets.VBox([angle_slider, output]))

    # Initialiser l'affichage
    update_angle({'new': angle_slider.value})

    return angle_slider

# Créer et afficher le widget
angle_widget = create_angle_widget()

# Fonction pour obtenir la valeur actuelle en radians
def get_angle_radians():
    return np.radians(angle_widget.value)
Taux de torsion (α) = 3.0020 rad/m
Rayon équivalent (R) = 0.2000 m
Rigidité en torsion (C) = 6.0358e+08 N·m²
# Calcul de la constante A pour le champ de torsion
A = 2 * C / (pi * R**4)  # A est un scalaire constant positif
print(f"Densité de contrainte de torsion (A_τ) = {A:.4e} N/m³")

# Définition des coordonnées spatiales
x = ufl.SpatialCoordinate(domain)

# Définition des composantes du champ de torsion (vecteur des Contraintes de Cauchy)
T1 = -A * x[1]  # Composante x du champ de torsion / Contrainte de cisaillement dans la direction x
T2 = A * x[0]   # Composante y du champ de torsion / Contrainte de cisaillement dans la direction y
T3 = ufl.as_ufl(0.0)  # Composante z nulle / Pas de contrainte dans la direction z (axe de torsion)

# Combinaison en un champ vectoriel de torsion
T = ufl.as_vector([T1, T2, T3])
print("Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)")
Densité de contrainte de torsion (A_τ) = 2.4016e+11 N/m³
Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)

Résolution#

Hide code cell content
# Définition des tenseurs de déformation et de contrainte
def epsilon(u):
    """Tenseur de déformation linéarisé"""
    return ufl.sym(ufl.grad(u))

def sigma(u):
    """Tenseur des contraintes (loi de Hooke)"""
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

# Définition des fonctions d'essai et de test
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Définition de la force volumique (nulle dans ce cas)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))

# Formulation variationnelle du problème
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx  # Forme bilinéaire

# Créez un widget Dropdown
surface_selector = widgets.Dropdown(
    options=[('Surface latérale Sl', 'Sl'), ('Surface inférieure S0', 'S0'), ('Surface supérieure Sh', 'Sh')],
    value='Sl',
    description='Surface:',
    disabled=False,
)

# Créez un widget de sortie pour afficher les messages
output = widgets.Output()

# Fonction pour mettre à jour la forme linéaire
def update_linear_form(change):
    selected_surface = change['new']
    # Mise à jour de la forme linéaire
    L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(surface_ids[selected_surface])
    
    # Mise à jour de l'affichage
    with output:
        output.clear_output(wait=True)
        print(f"Torsion appliquée sur {selected_surface}")

# Attacher la fonction de mise à jour au widget
surface_selector.observe(update_linear_form, names='value')

# Afficher le widget et la sortie
display(widgets.VBox([surface_selector, output]))

# Initialiser l'affichage
update_linear_form({'new': surface_selector.value})
# Créez un widget Dropdown pour la sélection de la condition aux limites
bc_selector = widgets.Dropdown(
    options=[
        ('Aucune', 'none'),
        ('Face inférieure S0', 'S0'),
        ('Face supérieure Sh', 'Sh'),
        ('Surface latérale Sl', 'Sl')
    ],
    value='S0',  # Valeur par défaut
    description='Encastrement:',
    disabled=False,
)

# Dictionnaire pour mapper les sélections aux conditions aux limites
bc_dict = {
    'none': [],
    'S0': [bc0_S0],
    'Sh': [bc0_Sh],
    'Sl': [bc0_Sl]
}

# Fonction pour mettre à jour le problème
def update_problem(change):
    selected_bc = change['new']
    problem = LinearProblem(a, L, bcs=bc_dict[selected_bc], 
                            petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
    
    # Mise à jour de l'affichage
    with output:
        output.clear_output(wait=True)
        print(f"Encastrement appliquée sur {selected_bc}")
            
# Attacher la fonction de mise à jour au widget
bc_selector.observe(update_problem, names='value')

# Afficher le widget
display(bc_selector, output)

# Initialisation du problème avec la valeur par défaut
update_problem({'new': bc_selector.value})
uh = problem.solve()

Visualisation des résultats#

Visualisation des déplacements#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Visualisation des rotations principales#

Visualisation des déplacements amplifiés#

Visualisation des déplacements normalisés et selon les axes principaux#

Visualisation des composantes du tenseur des déformations#

Visualisation des composantes du tenseur des contraintes#

Calcul des contraintes de von Mises#

s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))

V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)

Visualisation des contraintes de von Mises#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Analyses des résultats#

Calcul des valeurs maximales#

Hide code cell content
# Calcul du moment d'inertie polaire J pour une section circulaire
J = pi * (R ** 4) / 2

# Calcul de l'angle de torsion théorique
theta_theoretical = C * h / (mu * J)

# Extraction de la valeur maximale de la rotation RZ
max_rotation_z = np.max(np.abs(rotation_z.x.array))

# Extraction de la valeur maximale de la norme de rotation
max_rotation_norm = np.max(rotation_norm.x.array)

# Conversion en degrés pour une interprétation plus intuitive
max_rotation_z_degrees = np.degrees(max_rotation_z)

# Calcul de l'angle de torsion théorique
theoretical_twist = ang_radians * h  # ang est l'angle de torsion par unité de longueur, h est la hauteur du cylindre

# Comparaison pour RZ
relative_error_Z = abs(max_rotation_z - theoretical_twist) / theoretical_twist * 100

# Comparaison pour norme R
relative_error_N = (max_rotation_norm - theta_theoretical) / theta_theoretical * 100
Moment d'inertie polaire J : 2.513274e-03 m^4
Angle de torsion théorique θ : 3.001966 radians
Angle de torsion théorique θ : 172.00 degrés

Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 2.985258 radians
Rotation Z maximale simulée : 171.04 degrés
Erreur relative : 0.56%

Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 2.987141 radians
Norme de rotation maximale : 171.15 degrés
Erreur relative : -0.49%

Valeurs maximales des composantes de rotation :
RX max : 18.72 degrés
RX max  : 0.326659 radians
RY max : 18.33 degrés
RY max  : 0.320002 radians